{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Прогнозирование задержек вылетов рейсов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Введение"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наверное многие при отправлении в отпуск или командировку беспокоятся, как бы не опоздать на рейс. А если рейс утренний — то часто приходится искать компромисс между продолжительностью сна и ранним приездом в аэропорт. Бывает и так, что при планировании поездок приходится лететь с пересадками, и встаёт вопрос о том, хватит ли 40-50-60 минут, чтобы успеть на следующий рейс. Если бы при ответе на все эти вопросы у нас под рукой была информация, задержится ли вылет нашего рейса, скажем, на 15 минут — мы бы могли принять более взвешенное решение. Или просто могли бы спать на 15 минут дольше. В настоящем исследовании попробуем спрогнозировать задержку вылета рейса на 15 минут ради уменьшения количества синяков под глазами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Описание данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "За основу для анализа берём этот датасет с [данными](https://www.kaggle.com/usdot/flight-delays) о задержках внутренних рейсов американских авиакомпаний за 2015 год с Kaggle. К этому датасету добавим также [информацию](http://ourairports.com/data/) об аэропортах и взлётных полосах, а также историческую [информацию](https://www.ncdc.noaa.gov/isd/data-access) о погоде в аэропортах вылета и назначения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В исходном датасете, в файле ``flights.csv``, присутствуют следующие признаки:\n",
    "\n",
    "- Признаки, которые мы будем использовать:\n",
    "    - **YEAR**, **MONTH**, **DAY** — год, месяц и число, в которые выполнялся рейс\n",
    "    - **DAY_OF_WEEK** — день недели вылета\n",
    "    - **AIRLINE** — двухбуквенный код авиакомпании\n",
    "    - **ORIGIN_AIRPORT** — код аэропорта вылета\n",
    "    - **SCHEDULED_DEPARTURE** — запланированное время отбытия от гейта\n",
    "    - **DEPARTURE_DELAY** — задержка отбытия, ``DEPARTURE_TIME - SCHEDULED_DEPARTURE``\n",
    "    - **DISTANCE** — расстояние между аэропортами вылета и назанчения (в милях)\n",
    "- Признаки, которые мы не будем использовать для предсказаний:\n",
    "    - **FLIGHT_NUMBER** — номер рейса\n",
    "    - **TAIL_NUMBER** — бортовой номер самолёта, выполняющего рейс\n",
    "    - **DESTINATION_AIRPORT** — код аэропорта назначения\n",
    "    - **DEPARTURE_TIME** — реальное время отбытия от гейта\n",
    "    - **TAXI_OUT** — время, в течение которого самолёт осуществлял руление\n",
    "    - **WHEELS_OFF** — время вылета\n",
    "    - **SCHEDULED_TIME** — время полёта по плану\n",
    "    - **ELAPSED_TIME** — действительное время полёта (от гейта до гейта)\n",
    "    - **AIR_TIME** — время в воздухе (между взлётом и посадкой)\n",
    "    - **WHEELS_ON** — время посадки\n",
    "    - **TAXI_IN** — время руления в аэропорту назначения\n",
    "    - **SCHEDULED_ARRIVAL** — запланированное время прибытия к гейту в аэропорту назначения\n",
    "    - **ARRIVAL_TIME** — реальное время прибытия к гейту в аэропорту назначения\n",
    "    - **ARRIVAL_DELAY** — задержка прибытия, ``ARRIVAL_TIME - SCHEDULED_ARRIVAL``\n",
    "    - **DIVERTED** — самолёт приземлился не в аэропорту назначения\n",
    "    - **CANCELLED** — рейс был отменён\n",
    "    - **CANCELLATION_REASON** — код причины отмены рейса (A: авиакомпания, B: погода, C: National Airspace System, D: безопасность)\n",
    "    - **AIR_SYSTEM_DELAY** — величина задержки в минутах из-за National Airspace System\n",
    "    - **SECURITY_DELAY** — величина задержки в минутах в связи с вопросами безопасности\n",
    "    - **AIRLINE_DELAY** — величина задержки в минутах из-за авиакомпании\n",
    "    - **LATE_AIRCRAFT_DELAY** — величина задержки в минутах из-за задержки предыдущего рейса этого же самолёта\n",
    "    - **WEATHER_DELAY** — величина задержки в минутах из-за погоды\n",
    "\n",
    "Кроме того, в датасете содержатся ещё два файла — с информацией об аэропортах и расшифровкой кодов авиакомпаний. Информацию об аэропортах мы использовать не будем, поскольку есть альтернативный список аэропортов с большим количеством информации. Файл ``airlines.csv`` имеет следующие поля:\n",
    "\n",
    "- **IATA_CODE** — код авиакомпании (соответствует полю ``AIRLINE`` из файла со списком рейсов)\n",
    "- **AIRLINE** — название авиакомпании"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Первичный анализ данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import datetime as dt\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загружаем датасет, смотрим типы данных:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('Data/US2015/flights.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, всего в датасете есть информация о 5.819kk совершённых рейсов."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выкинем лишние поля"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop(['FLIGHT_NUMBER', 'TAIL_NUMBER', 'DEPARTURE_TIME', 'TAXI_OUT', 'WHEELS_OFF', 'SCHEDULED_TIME', \n",
    "         'ELAPSED_TIME', 'AIR_TIME', 'WHEELS_ON', 'TAXI_IN', 'SCHEDULED_ARRIVAL', 'ARRIVAL_TIME',\n",
    "         'ARRIVAL_DELAY', 'DIVERTED', 'CANCELLED', 'CANCELLATION_REASON'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Целевая переменная"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сформируем целевую переменную из поля ``DEPARTURE_DELAY``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['target'] = df['DEPARTURE_DELAY'] >= 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('{:.2%} рейсов было задержано на 15 и более минут'.format(df['target'].sum() / len(df)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Поскольку целевая переменная может принимать 2 значения, будем решать задачу классификации"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Основные причины задержек"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим распределение причин задержек рейсов (долю рейсов, в которых данная причина указана)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(df.loc[df['target'], \n",
    "        ['AIR_SYSTEM_DELAY', 'SECURITY_DELAY', 'AIRLINE_DELAY',\n",
    "         'LATE_AIRCRAFT_DELAY', 'WEATHER_DELAY']] > 0).sum() / df['target'].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Во-первых, видим, что сумма всех причин задержек больше 1, то есть одновременно несколько факторов могут приводить к задержке. Во-вторых, наиболее вероятными причинами задержки является задержка прибытия предыдущего рейса или что-то связанное с авиакомпанией (подготовка самолёта, ожидание пасажиров или экипажа, удаление с борта буйных пасспажиров и т.п., на [сайте](http://aspmhelp.faa.gov/index.php/Types_of_Delay) Федерального Авиационного Агентства есть список с описанием причин)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Аэропорты и авиакомпании"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, какие аэропорты имеются в датасете:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_airports = set(df['ORIGIN_AIRPORT'].unique().tolist() + df['DESTINATION_AIRPORT'].unique().tolist())\n",
    "print('В датасете присутствует информация о {} аэропортах:'.format(len(unique_airports)))\n",
    "print(unique_airports)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Цифровые коды аэропортов — это какая-то проблема с данными, другой вид кодировки. Находим следующий [список](https://github.com/srcole/flightdelay/tree/master/airportcodes) аэропортов с соответствием трёхбуквенных и цифровых кодов. По ссылке два списка аэропортов — с полями **Code** (в одном файле это IATA код, в другом цифровой) и **Description** (название аэропорта в обоих файлах). Заменим цифровые коды на буквенные IATA коды."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "codes_num = pd.read_csv('Data/raw/L_AIRPORT_ID.csv')\n",
    "codes_iata = pd.read_csv('Data/raw/L_AIRPORT.csv')\n",
    "\n",
    "class Vocab(dict):\n",
    "    def __missing__(self, key):\n",
    "        return key\n",
    "\n",
    "fix_vocab = Vocab(\n",
    "    pd.merge(codes_num, codes_iata, on='Description')[['Code_x', 'Code_y']] \\\n",
    "            .astype(str).set_index('Code_x').to_dict()['Code_y'])\n",
    "\n",
    "df['ORIGIN_AIRPORT'] = df['ORIGIN_AIRPORT'].astype(str).map(fix_vocab)\n",
    "df['DESTINATION_AIRPORT'] = df['DESTINATION_AIRPORT'].astype(str).map(fix_vocab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Попробуем ещё раз:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_airports = set(df['ORIGIN_AIRPORT'].unique().tolist() + df['DESTINATION_AIRPORT'].unique().tolist())\n",
    "print('В датасете присутствует информация о {} аэропортах:'.format(len(unique_airports)))\n",
    "print(unique_airports)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Гораздо лучше. Теперь посмотрим на авиакомпании."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Авиакомпании, представленные в датасете: {} (всего {} штук)'.format(df['AIRLINE'].unique().tolist(),\n",
    "                                                                           len(df['AIRLINE'].unique().tolist())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Здесь всё правильно, сокращения совпадают с приведёнными в файле ``airlines.csv``. Не будем доставать оттуда полные названия компаний, чтобы потом их было проще кодировать для обучения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Визуальный анализ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Распределение вероятностей задержек по времени"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для удобства преобразуем данные, связанные с датой и временем к ``datetime`` формату."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_time(x):\n",
    "    hour = int(x[:2])\n",
    "    minute = int(x[-2:])\n",
    "    return dt.time(hour, minute)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df['SCHEDULED_DEPARTURE'] = df['SCHEDULED_DEPARTURE'].astype(str).str.zfill(4).apply(get_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df['DATE'] = pd.to_datetime(df[['YEAR', 'MONTH', 'DAY']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%time\n",
    "# df['DATETIME'] = df.apply(lambda x: dt.datetime.combine(x['DATE'], x['SCHEDULED_DEPARTURE']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby('DAY_OF_WEEK')['target'].agg(['sum', 'count'])\n",
    "xs = gr.index.values\n",
    "ys = gr['sum'] / gr['count']\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "sns.barplot(xs, ys);\n",
    "\n",
    "plt.xlabel('День недели');\n",
    "plt.ylabel('Вероятность задержки рейса');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Существенной зависимости от дня недели нет."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby(pd.Grouper(key='DATE', freq='D'))['target'].agg(['sum', 'count'])\n",
    "ys = gr['sum'] / gr['count']\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(16, 6))\n",
    "plt.plot(np.arange(len(ys)), ys);\n",
    "plt.xlabel('День года');\n",
    "plt.ylabel('Вероятность задержки рейса');\n",
    "plt.tick_params('y', colors='C0')\n",
    "plt.ylim(0, 0.45)\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ys = gr['count']\n",
    "plt.plot(np.arange(len(ys)), ys, lw=3, color='C1');\n",
    "plt.ylim(0, 20000);\n",
    "plt.ylabel('Число рейсов');\n",
    "plt.tick_params('y', colors='C1');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Интересно, что в районе нового года и рождества вероятность задержки рейса существенно растёт, хотя количество выполняемых рейсов остаётся на обычном уровне. Возможно это связано с тем, что в аэропортах остаётся работать очень небольшое количество работников."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "То же самое, но по неделям:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby(pd.Grouper(key='DATE', freq='W'))['target'].agg(['sum', 'count'])\n",
    "ys = gr['sum'] / gr['count']\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(16, 6))\n",
    "plt.bar(np.arange(len(ys)), ys);\n",
    "plt.xlabel('День года');\n",
    "plt.ylabel('Вероятность задержки рейса');\n",
    "plt.tick_params('y', colors='C0')\n",
    "plt.ylim(0, 0.45)\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ys = gr['count']\n",
    "plt.plot(np.arange(len(ys)), ys, lw=3, color='C1');\n",
    "plt.ylim(0, 120000);\n",
    "plt.ylabel('Число рейсов');\n",
    "plt.tick_params('y', colors='C1');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Введём дополнительно время вылета в часах (дробных и целых):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['DEP_HOUR_FLOAT'] = df['SCHEDULED_DEPARTURE'].apply(lambda x: x.hour + x.minute/60)\n",
    "df['DEP_HOUR_INT'] = df['DEP_HOUR_FLOAT'].astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby('DEP_HOUR_INT')['target'].agg(['sum', 'count'])\n",
    "xs = gr.index.values\n",
    "ys = gr['sum']/gr['count']\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.bar(xs, ys);\n",
    "\n",
    "plt.xlabel('Час дня');\n",
    "plt.ylabel('Вероятность задержки');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что в течение дня вероятность задержки рейса существенно меняется. Возможно это связано с тем, что ближе к концу дня задержки рейсов накапливаются и задержка предыщих рейсов приводит к задержке последующих. В таком случае количество вероятность задержки данного рейса будет тем больше, чем больше рейсов уже вылетело из данного аэропорта на текущий момент. Посмотрим:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df['FLIGHTS_BEFORE'] = df.groupby(['DATE', 'ORIGIN_AIRPORT'])['DATE'].transform(lambda x: list(range(len(x))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим вероятность задержки по логарифмически распределённым корзинам:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_id = np.digitize(df['FLIGHTS_BEFORE'],\n",
    "                     bins=np.logspace(0, 10, 11, base=2))\n",
    "\n",
    "df['FLIGHT_BEFORE_BIN'] = bin_id\n",
    "\n",
    "for b in np.unique(bin_id):\n",
    "    flt = bin_id == b\n",
    "    plt.bar(b, np.sum(df.loc[flt, 'target']) / np.sum(flt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Действительно — чем больше рейсов вылетело из аэропорта, тем выше вероятность задержки текущего рейса."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построем распределение вероятности задержки одновременно от часа вылета и корзины по количествую предшествующих вылетов, в которую попадает данный рейс."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_meshgrid(x, y, z):\n",
    "    ux = np.unique(x)\n",
    "    uy = np.unique(y)\n",
    "    [X, Y] = np.meshgrid(ux, uy)\n",
    "    \n",
    "    def getz(x1, y1):\n",
    "        fltx = x == x1\n",
    "        flty = y == y1\n",
    "        try:\n",
    "            idx = (fltx & flty).tolist().index(True)\n",
    "        except ValueError:\n",
    "            idx = 0\n",
    "        return idx\n",
    "    \n",
    "    Z = np.array([z[getz(x1, y1)] for x1, y1 in zip(X.flatten(), Y.flatten())])\n",
    "    Z.shape = X.shape\n",
    "    \n",
    "    return (X, Y, Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby(['DEP_HOUR_INT', 'FLIGHT_BEFORE_BIN'])['target'].agg(['sum', 'count']).reset_index()\n",
    "xs = gr['DEP_HOUR_INT']\n",
    "ys = gr['FLIGHT_BEFORE_BIN']\n",
    "zs = gr['sum'] / gr['count']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, Y, Z = gen_meshgrid(xs, ys, zs)\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.pcolor(X, Y, Z, cmap='magma')\n",
    "plt.colorbar();\n",
    "plt.clim(0, 0.3);\n",
    "plt.xlabel('Час дня');\n",
    "plt.ylabel('Количество предшествующих вылетов (log2)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что у тех вылетающих вечером рейсов, перед которыми было мало других вылетов, ниже вероятность задержки."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим также как меняется вероятность задержки в зависимости от дня недели и времени суток."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby(['DEP_HOUR_INT', 'DAY_OF_WEEK'])['target'].agg(['sum', 'count']).reset_index()\n",
    "xs = gr['DEP_HOUR_INT']\n",
    "ys = gr['DAY_OF_WEEK']\n",
    "zs = gr['sum'] / gr['count']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, Y, Z = gen_meshgrid(xs, ys, zs)\n",
    "\n",
    "plt.figure(figsize=(16, 4))\n",
    "sns.heatmap(Z, cmap='magma');\n",
    "plt.xlabel('Час дня');\n",
    "plt.ylabel('День недели');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что в целом вероятности довольно однородны в течение недели. Немного выделяются вечер четверга (вероятность задержки выше) и вечер субботы (вероятность ниже)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Распределение задержек по авиакомпаниям"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на проценты задержек у разных авикомпаний."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gr = df.groupby('AIRLINE')['target'].agg(['sum', 'count'])\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "xs = gr.index.values\n",
    "ys = gr['sum'] / gr['count']\n",
    "sns.barplot(xs, ys);\n",
    "\n",
    "plt.xlabel('Код авиакомпании');\n",
    "plt.ylabel('Вероятность задержки');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Влияние дальности перелёта"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.histogram(df['DISTANCE'], bins=40)[1]\n",
    "bin_id = np.digitize(df['DISTANCE'], bins=bins)\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(12, 6))\n",
    "for bid, b in zip(np.unique(bin_id), bins):\n",
    "    flt = bin_id == bid\n",
    "    plt.bar(b, np.sum(df.loc[flt, 'target']) / np.sum(flt), width=50, color='C0')\n",
    "\n",
    "plt.ylabel('Вероятность задержки');    \n",
    "\n",
    "ax1.twinx()\n",
    "for bid, b in zip(np.unique(bin_id), bins):\n",
    "    flt = bin_id == bid\n",
    "    plt.bar(b+50, np.sum(flt), width=50, color='C1')\n",
    "    \n",
    "plt.yscale('log')    \n",
    "\n",
    "plt.xlabel('Дальность перелёта');\n",
    "plt.ylabel('Количество перелётов в корзине');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Интересно, что для рейсов дальностью больше 3000 миль, резко увеличивается вероятность задержек. По-видимому такие перелёты всегда сопровождаются более тщательной проверкой и подготовкой самолёта. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(df.loc[(df['DISTANCE'] > 3000) & (df['target']),\n",
    "        ['AIR_SYSTEM_DELAY', 'SECURITY_DELAY', 'AIRLINE_DELAY',\n",
    "         'LATE_AIRCRAFT_DELAY', 'WEATHER_DELAY']] > 0).sum() / ((df['DISTANCE'] > 3000) & (df['target'])).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И действительно, в причинах задержек дальних рейсов единолично доминируют причины, связанные с авиакомпанией."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Подготовка фич "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Достанем ещё несколько фич из таблицы с аэропортами. Добавим координаты аэропорта, его размер, высоту над уровнем моря, наличие планового обслуживания и регион, в котором он находится (страна и штат США)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airport_df = pd.read_csv('Data/raw/airports.csv').set_index('id')\n",
    "\n",
    "lat_vocab, lon_vocab, type_vocab = {}, {}, {}\n",
    "elev_vocab, region_vocab, service_vocab = {}, {}, {}\n",
    "for airport in unique_airports:\n",
    "    flt = (airport_df['iata_code'] == airport) & (airport_df['type'] != 'closed')\n",
    "    lat_vocab[airport] = airport_df.loc[flt, 'latitude_deg'].values[0]\n",
    "    lon_vocab[airport] = airport_df.loc[flt, 'longitude_deg'].values[0]\n",
    "    type_vocab[airport] = airport_df.loc[flt, 'type'].values[0]\n",
    "    elev_vocab[airport] = airport_df.loc[flt, 'elevation_ft'].values[0]\n",
    "    region_vocab[airport] = airport_df.loc[flt, 'iso_region'].values[0]\n",
    "    service_vocab[airport] = airport_df.loc[flt, 'scheduled_service'].values[0]\n",
    "\n",
    "df['ORIGIN_LAT'] = df['ORIGIN_AIRPORT'].map(lat_vocab)\n",
    "df['ORIGIN_LON'] = df['ORIGIN_AIRPORT'].map(lon_vocab)\n",
    "df['ORIGIN_SIZE'] = df['ORIGIN_AIRPORT'].map(type_vocab)\n",
    "df['ORIGIN_ELEV'] = df['ORIGIN_AIRPORT'].map(elev_vocab)\n",
    "df['ORIGIN_REGION'] = df['ORIGIN_AIRPORT'].map(region_vocab)\n",
    "df['ORIGIN_SERVICE'] = df['ORIGIN_AIRPORT'].map(service_vocab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Добавим в датафрейм день года."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['DAY_OF_YEAR'] = pd.to_datetime(df['DATE']).dt.dayofyear"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Выбор модели и метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Поскольку решается задача бинарной классификации с несбалансированными классами, будем пользоваться метрикой ROC AUC. Обучать будем модели на основе деревьев — ``XGBoostClassifier``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "from sklearn.base import BaseEstimator, TransformerMixin\n",
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "from sklearn.pipeline import Pipeline, FeatureUnion\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from xgboost import XGBClassifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для начала соберём пайплайн для предобработки данных и обучения модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ItemSelector(BaseEstimator):\n",
    "    def __init__(self, colname):\n",
    "        self.colname = colname\n",
    "        pass\n",
    "    \n",
    "    def fit(self, x, y=None):\n",
    "        return self\n",
    "    \n",
    "    def transform(self, x, y=None):\n",
    "        if x[self.colname].dtype == 'object':\n",
    "            return x[self.colname].values\n",
    "        else:\n",
    "            return x[self.colname].values.astype('float64').reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_pipeline = Pipeline([\n",
    "    ('union', FeatureUnion([\n",
    "        ('dep_hour', ItemSelector('DEP_HOUR_FLOAT')),\n",
    "        \n",
    "        ('orig_airport', Pipeline([\n",
    "            ('selector', ItemSelector('ORIGIN_AIRPORT')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ])),\n",
    "        \n",
    "        ('airline', Pipeline([\n",
    "            ('selector', ItemSelector('AIRLINE')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ])),\n",
    "        \n",
    "        ('distance', ItemSelector('DISTANCE')),\n",
    "    ])),\n",
    "    \n",
    "    ('xgb', XGBClassifier(random_state=241))\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``ItemSelector`` выделяет из всего датафрейма нужную фичу и отдаёт её для дальнейшей обработки, либо просто для добавления в конечную матрицу признаков-наблюдений.\n",
    "\n",
    "``CountVectorizer`` в данном случае работает как OneHotEncoder для кода аэропорта, кода авиакомпании, размера аэропорта, наличия сервиса в аэропорту. Для региона (формат *КОДСТРАНЫ-ШТАТ*) ``CountVectorizer`` достаёт отдельно страну и отдельно штат."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Разделение на обучающую и тестовую выборки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вообще говоря, мы бы хотели научиться предсказывать будущие задержки рейсов. Однако, поскольку данные у нас есть только за один год, в течение которого наблюдаются сезонные колебания, разделять выборку на обучающую и тестовую по времени будет неверно — в обучающей выборке не будет особенностей, связанных, скажем, с рождественскими праздниками. Можно разделять по аэропортам (например, разделив их на западную и восточную части), но в таком случае также из обучающей выборки потеряются особенности, представленные в тестовой. Поэтому будем разобъём весь датасет на группы с уникальными парами (день, аэропорт) и разделим датасет по ним. В таком случае в обучающей выборке у нас будут присутствовать все даты и все аэропорты, но в то же время, мы избежим утечек данных о конкретном аэропорте в конкретный день."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_groups = list(df.groupby(['DAY_OF_YEAR', 'ORIGIN_AIRPORT']).groups.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "np.random.seed(27)\n",
    "train_size = int(len(unique_groups) * 0.4)\n",
    "train_groups = random.sample(unique_groups, train_size)\n",
    "test_groups = list(set(unique_groups) - set(train_groups))\n",
    "df_train = df.set_index(['DAY_OF_YEAR', 'ORIGIN_AIRPORT'], drop=False).loc[train_groups].reset_index(drop=True)\n",
    "df_test = df.set_index(['DAY_OF_YEAR', 'ORIGIN_AIRPORT'], drop=False).loc[test_groups].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Доля задержек в обучающей выборке: {:.2%}'.format(df_train['target'].sum() / len(df_train)))\n",
    "print('Доля задержек в тестовой выборке: {:.2%}'.format(df_test['target'].sum() / len(df_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Baseline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обучим базовую модель"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "xgb_pipeline = xgb_pipeline.fit(df_train, df_train['target']);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = xgb_pipeline.predict_proba(df_test)[:,1]\n",
    "print('ROC AUC базовой модели: {:.4f}'.format(roc_auc_score(df_test['target'], y_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Добавим дополнительные признаки (количество вылетов до текущего рейса, признаки, связанные с аэропортом)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_pipeline = Pipeline([\n",
    "    ('union', FeatureUnion([\n",
    "        ('dep_hour', ItemSelector('DEP_HOUR_FLOAT')),\n",
    "        \n",
    "        ('orig_airport', Pipeline([\n",
    "            ('selector', ItemSelector('ORIGIN_AIRPORT')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ])),\n",
    "        \n",
    "        ('airline', Pipeline([\n",
    "            ('selector', ItemSelector('AIRLINE')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ])),\n",
    "        \n",
    "        ('distance', ItemSelector('DISTANCE')),\n",
    "        \n",
    "        ('flights_before', ItemSelector('FLIGHTS_BEFORE')),\n",
    "        \n",
    "        ('dayofyear', ItemSelector('DAY_OF_YEAR')),\n",
    "        \n",
    "        ('dayofweek', ItemSelector('DAY_OF_WEEK')),\n",
    "        \n",
    "        ('airport_size', Pipeline([\n",
    "            ('selector', ItemSelector('ORIGIN_SIZE')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ])),\n",
    "        \n",
    "        ('lat', ItemSelector('ORIGIN_LAT')),\n",
    "        \n",
    "        ('lon', ItemSelector('ORIGIN_LON')),\n",
    "        \n",
    "        ('elevation', ItemSelector('ORIGIN_ELEV')),\n",
    "        \n",
    "        ('service', Pipeline([\n",
    "            ('selector', ItemSelector('ORIGIN_SERVICE')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ])),\n",
    "        \n",
    "        ('region', Pipeline([\n",
    "            ('selector', ItemSelector('ORIGIN_REGION')),\n",
    "            ('ohe', CountVectorizer())\n",
    "        ]))\n",
    "    ])),\n",
    "    \n",
    "    ('xgb', XGBClassifier(random_state=241, n_estimators=100))\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Кросс-валидация и настройка гиперпараметров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проведём кросс-валидацию на 2 фолдах. Поскольку у нас очень много наблюдений в обучащей выборке, мы и без использования StratifiedKFold должны получить примерно одинаковый % задержек на разных фолдах. Будем искать лучшее значение параметра ``learning_rate``, который отвечает за то, какой вклад в общую модель даёт каждое из деревьев."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {'xgb__learning_rate': [0.001, 0.01, 0.1]}\n",
    "gs = GridSearchCV(xgb_pipeline, params, scoring='roc_auc', cv=2, verbose=1, return_train_score=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "gs = gs.fit(df_train, df_train['target'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('ROC AUC лучше модели на кросс-валидации: {:.4f}'.format(gs.best_score_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим кривые валидации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(params['xgb__learning_rate'], gs.cv_results_['mean_test_score']);\n",
    "plt.fill_between(params['xgb__learning_rate'],\n",
    "                gs.cv_results_['mean_test_score'] - gs.cv_results_['std_test_score'],\n",
    "                gs.cv_results_['mean_test_score'] + gs.cv_results_['std_test_score']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Конечно надо изучать больше параметров, больше значений для каждого и на большем количестве фолдов, но тогда можно не успеть к дедлайну."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, что получилось на отложенной выборке."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = gs.predict_proba(df_test)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('ROC AUC базовой модели: {:.4f}'.format(roc_auc_score(df_test['target'], y_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Значение метрики на отложенной выборке очень близко к результату кросс-валидации. По сравнению с бейзлайном результат улучшился, но фактически только за счёт новых фич, поскольку на кросс-валидации лучшими оказались дефолтные параметры."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Выводы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть существенное пространство для улучшения результата по многим направлениям. Можно добавить фичу, которая будет содержать количество вылетов за небольшой период перед текущим рейсом, что будет отражать загруженность аэропорта в текущий момент. Есть открытые [данные](https://www.ncdc.noaa.gov/isd/data-access) о погоде вблизи аэропортов с часовой периодичностью — эту информацию можно использовать, чтобы оценивать вероятность задержки из-за плохой погоды. Безусловно необходимо провести более тщательную настройку гиперпараметров — в первую очередь количества деревьев в модели (``n_estimators``), максимальной глубины деревьев (``max_depth``)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
